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Abstract. We investigate the long time stability in Nekhoroshev's sense for the Sun- Jupiter— 
^ Saturn problem in the framework of the problem of three bodies. Using computer algebra in 

Xf^ . order to perform huge perturbation expansions we show that the stability for a time comparable 

CTN \ with the age of the universe is actually reached, but with some strong truncations on the 

perturbation expansion of the Hamiltonian at some stage. An improvement of such results is 
currently under investigation. 



1. Introduction 

The stability of the Solar System is a classical, long standing and challenging problem, 
already pointed out by Newton. In this article we revisit the problem in the light of the 
theorems of Kolmogorov and of Nekhoroshev, with the aim of proving that they apply 
to the problem of three bodies with the masses and orbital parameters of Jupiter and 
Saturn. 

Let us briefly recall the historical development of our knowledge. After 1954 a 
possible solution was suggested by the celebrated theorem of Kolmogorov I stating 
the existence of a large measure set of invariant tori for a nearly integrable Hamiltonian 
system, e.g., the planetary system when the mutual perturbation of the planet is taken 
into account. The relevance of Kolmogorov's result for the planetary problem has been 
soon emphasized by Arnoldt^^ and Moserl"^"^!. In particular Arnold worked out a proof 
taking into account the degeneration of the unperturbed Hamiltonian which occurs in 
the planetary case^^^. On the other hand, Moser first gave a proof for the case of an area 
preserving mapping of an annulus^'^^] , and a few years later pointed out that the theorem 
of Kolmogorov implies that the classical Lindstedt series are actually convergent ["^^1. 
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As a matter of fact it was soon remarked by Henon that the appHcation of the Kol- 
mogorov's theorem to the planetary motions is not straightforward, due to the condition 
that the masses of the planets should be small enough. Indeed, the available estimates 
could only assure the applicability, e.g., to the problem of three bodies, if the masses of 
the planets are less than that of a proton. On the other hand, numerical integrations of 
the full Solar System over a time span of billions of years have shown that the orbits of 
the inner planets exhibit a chaotic evolution which is incompatible with the quasi pe- 
riodic motion predicted by Kolmogorov's theoremt^lt^^lt^-^^. Furthermore the subsystem 
of the major planets (i.e. Jupiter, Saturn, Uranus and Neptune) shows a very small pos- 
itive Lyapunov exponent t^'^l, once again, this cannot fit with a motion on an invariant 
torus. 

A second approach was suggested by Moscrl*^^] and Littlewood^^'^lt^^l and fully stated 
by Nekhoroshevt"^^! ["^^1 with his celebrated theorem on exponential stability. According 
to this theorem the time evolution of the actions of the system (which in the planetary 
case are actually related to the semimajor axes, the inclinations and the eccentricities) 
remains bounded for a time exponentially increasing with the inverse of the pertur- 
bation parameter. Thus, although the possibility of a chaotic motion is not excluded, 
nevertheless a dramatic change of the orbits should not occur for such a long time, and 
it may be conjectured that such a time exceeds the age of the Solar System itself. But 
also in this case the problem of the applicability of the theorem still persists, since the 
analytical estimates based on Nekhoroshev's formulation or other analytical proofs give 
ridiculous estimates for the size of the masses of the planets. 

In recent years the estimates for the applicability of both Kolmogorov's and 

Nekhoroshev's theorems to realistic models of some part of the Solar System have been 
improved by some authors. For example the applicability of Nekhoroshev's theorem 
to the stability of the Trojan asteroids in the vicinity of the triangular Lagrangian 
points has been investigated by Giorgilli et al.I^^'Wt-^^l), Efthimiopoulos et al.t^l and by 
Lhotka et al.t^^', the connection between Nekhoroshev's theorem and Arnold diffusion 
has been considered by Efthimiopoulost^*^!; the applicability of KAM theorem has been 
studied by Robutelt"^^!, Fejozl^^l, Celletti et alJ^\ Gabern et al.t^^l and by Locatelli and 
GiorgilliP^lP^]. In particular in the latter two articles the Sun- Jupiter-Saturn (hereafter 
SJS) system is investigated, and evidence is produced that an invariant torus exists in 
the vicinity of the initial data of Jupiter and Saturn, at least in the approximation of 
the general problem of three bodies. 

In the present article we study the stability in Nekhoroshev's sense of the neigh- 
bourhood of the invariant torus for the SJS system. The aim is to give evidence, with 
help of a computer-assisted calculation, that the size of the neighbourhood of the invari- 
ant torus for which exponential stability holds for a time interval as long as the age of 
the universe is big enough to contain the actual initial data of Jupiter and Saturn. We 
should say that such an ambitious goal is still out of our actual possibilities. However, 
we show that our methods should allow us to achieve our program provided a sufficient 
computer power will be available in the next future and a further refinement of our 
approximation methods will be worked out. This is work for the next future. 
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2. Theoretical framework 



The basis of our approach is the investigation of the stabihty of a neighbourhood of 
an invariant Kolmogorov's torus. To this end let us briefly recall the statement of Kol- 
mogorov's theorem 

Theorem 1: Consider a canonical system with Hamiltonian 
(1) Hip, q) = h{p) + £/(p, q) . 

Let us assume that the unperturbed part of the Hamiltonian is non-degenerate, i.e., 
det (^ Qp.Qp^ ^ 7^ 0, and that p* e R"' is such that the corresponding frequencies u = 
^{p*) satisfy a Diophantine condition, i.e., 

\{k,uj)\>^\k\-^ VO^/ceZ'^ 

with some constants 7 > and t > n — 1 . Then for e small enough the Hamiltonian (1 ) 
admits an invariant torus carrying quasiperiodic motions with frequencies u). The invari- 
ant torus lies in a e-neighbourhood of the unperturbed torus {{p,q) '■ P = P*, 9 G T"^}. 

The question is about the dynamics in the neighbourhood of the invariant torus. 
In order to discuss this point we need a few technical details about the Kolmogorov's 
proof method. The key points, clearly outlined in the original short note [19], are the 
following. First, one picks an unperturbed invariant torus p* for the Hamiltonian (1) 
characterized by diophantine frequencies a;, and expands the Hamiltonian in power series 
of the actions p in the neighbourhood of p* . Thus (with a translation moving p* to the 
origin of the actions space) one gives the initial Hamiltonian the form 



(2) H{p, q) = {io,p) + sA{q) + e{B{q),p) + ^{Cp,p) + 0(p^) 

where C = g^-^ip*) is a symmetric matrix, and A{q) and (^B{q),p) are the terms 

independent of p and linear in p in the power expansion of the perturbation f{p, q) , 
respectively. The quadratic part in OijP') is of order £, too. The next step consists in 
performing a near the identity canonical transformation which gives the Hamiltonian 
the Kolmogorov's normal form 

(3) H'{p',q') = {uj,p')+0{p'^) . 

As Kolmogorov points out, the invariance of the torus = is evident, due to the par- 
ticular form of the normalized Hamiltonian. The whole process requires a composition 
of an infinite sequence of transformations, and the most difficult part is to prove the con- 
vergence of such a sequence. The point which is of interest to us is that the transformed 
Hamiltonian (3) is analytic in a neighbourhood of the invariant torus p' = . 

Let us emphasize that the analytical form of the Hamiltonian (3) is quite similar 
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to that of a Hamiltonian in the neighbourhood of an eUiptic equihbrium, namely 

1 

where the dots stand for terms of degree larger than 2 in the Taylor expansion. For, 
introducing the action-angle variables p, q via the usual canonical transformation Xj = 
^2pj cos Qj , Uj = -s/^pj sin qj , the latter Hamiltonian takes essentially the form (3). 
Thus the exponential stability of the invariant torus p' = may be proved using the 
theoretical scheme that works fine in the case of an elliptic equilibrium, e.g., in the case 
the triangular Lagrangian points. 

As a matter of fact, a much stronger result holds true, namely that the invariant 
torus is superexponentially stable, as stated in [29] and [15]. However, a computer 
assisted method for the theory of superexponential stability seems not be currently 
available, so we limit our study to the exponential stability in Nekhoroshev's sense. 



3. Technical tools 

Let us now come to the improvement of the estimates for the applicability of the theo- 
rems of Kolmorogov and Nekhoroshev. The key point is to use an explicit construction 
of the normal form up to a finite order with algebraic manipulation in order to reduce 
the size of the perturbation, and then apply a suitable formulation of the theorems. 

Let us explain this point by making reference to the theorem of Kolniogorov. Start- 
ing with the Hamiltonian (2), we perform a finite number, r say, of normalization steps 
in order to give the Hamiltonian the normal form up to order r 

(4) H^'^\p,q) = {0J,p) + ^{Cp,p)+e^A^^\q)+e^{B'^^\q),p)+n'^^\p,q) 

with TZ^'^\p, q) = Cdpp), so that the perturbation is now of order e^. 

To this end we implement the normalization algorithm for the normal form of 
Kolmogorov step by step in powers of e, as in the traditional expansions in Celestial 
Mechanics. The full justification of such a procedure, including the convergence proof, 
is given, e.g., in [16] and [17]. The resulting Hamiltonian has still the form (2), with, 
however, e replaced by e'^ . Thus, a straightforward application of the theorem reads, in 
rough terms: if < £*, then an invariant torus exists. The power r may considerably 
improve the estimate of the threshold for the applicability of the theorem. This approach 
has been translated in a computer assisted rigorous proof, which has been successfully 
applied to a few simple modelst^lt^^^^^^^). 

Let us now come to the part concerning the estimate of the stability time which 
is the main contribution of the present note. To this end we remove from the Hamilto- 
nian (4) all the contributions which are independent of or linear in the actions p, namely 
the terms e which are small, thus obtaining a reduced Hamil- 

tonian in Kolmogorov's normal form. Moreover, we expand the perturbation TiS^'^ {jp, q) 
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in power series of p and Fourier series of q, thus getting a Hamiltonian in the form 

(5) H{p, q) = {uj,p) + Hi{p, q) + H2{p, q) + . . . • 

where Hg (p, q) is a homogeneous polynomial of degree s + 1 in the actions p and a 
trigonometric series in the angles q. Here, the upper index r of H has been removed 
because it is now meaningless, since we use the latter Hamiltonian as an approximation 
of the Kolmogorov's normal form. 

On this Hamiltonian we perform a Birkhoff normalization up to a finite order, 
that we denote again by r although it has no relation with the order of Kolmogorov's 
normalization used above. Thus we get a Birkhoff normalized Hamiltonian 

H = {uj,p) + Zxip) + . . . + Zr{p) + Trip, q), 

with J>(p, a power series in p starting with terms of degree r + 2. We omit the 
details about this part of the calculation, since there are a number of well known formal 
algorithms that do the job. We concentrate instead on the quantitative estimates. 

Let us introduce a norm for a function /(p, q) = X]|z|=s fceZ" /i.fcP^e'^'^'^^ which is a 
homogeneous polynomial of degree s in the actions p. Precisely define 

(6) 11/11= E 

\i\=s , feez" 

Moreover consider the domain 

(7) A, = {peR'', \pj\<g, j = l, ... ,n} . 
Then we have 

|/(p,g)|< 11/11^^ forpe A,, geT- . 

Let now p{0) G A^^ with < Q- Then we have p{t) G A^ for \t\ < T, where T is the 
escape time from the domain A^. This is the quantity that we want to evaluate. To this 
end we use the elementary estimate 

(8) \pit) -p(0)| < \t\ ■ sup \p\ < \t\ ■ \\{p,T}\\g'-+^ . 

\p\<e 

The latter formula allows us to find a lower bound for the escape time from the domain 
Ag, namely 

(9) r{QQ,Q,r)^ 



P,F}\\q-+^ 

which however depends on Q and r. We emphasize that in a practical application, 
e.g., to the SJS system, is fixed by the initial data, while g and r are left arbitrary. 
Thus we try to find an estimate of the escape time T{qq) depending only on the physical 
parameter ^o- To this end we optimize t(^o, Qi f) with respect to q and r, proceeding as 
follows. First we keep r fixed, and remark that the function t{qq^ g, r) has a maximum 
for 

r + 2 
r + 1 
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This gives an optimal value of ^ as a function of and r, and so a new function 

riQo,r) = sup T{go,Q,r) 

Q>QO 

which is actually computed by putting the optimal value g = ^o(^ + 2)/(r + 1) in the 
expression above for t^qq^q^t). Next we look for the optimal value ropt of r, which 
maximizes f^Qo^r) when r is allowed to change. That is, we look for the quantity 

T{qo) = maxf(^o,r) , 

r>l 

which is our best estimate of the escape time, depending only on the initial data. We 
define the latter quantity as the estimated stability time. Here we need some further 
considerations in order to convince the reader that the maximum in the r.h.s. of the latter 
formula actually exists. This follows from the asymptotic properties of the Birkhoff's 
normal form. For, according to the available analytical estimates based on Diophantine 
inequalities for the frequencies, the norm || {p, J>} || in the denominator of (9) is expected 
to grow as (r!)"^, n being the number of degrees of freedom. Thus, for small enough 
the denominator ^r}||^'o reaches a minimum for some r"^ ~ ^/ Qo: which means that 
the wanted maximum actually exists, thus providing the optimal value ropt • We also 
remark that although no proof exists that the analytical estimates are optimal, accurate 
numerical investigations based on explicit expansions show that the r! growth of the 
norms actually shows up (see, e.g., [7] and [8]). Working out an analytical evaluation of 
the stability time on the basis of these considerations leads to an exponential estimate 
of Nekhoroshev's type for T{go) (see, e.g., [13]). Here we replace the analytical estimates 
with an explicit numerical optimization of f{go,r) by just calculating it for increasing 
values of r until the maximum is reached. 

Our aim is to perform the procedure above by using computer algebra. Thus some 
truncation of the functions must be introduced in order to implement the actual cal- 
culation. The most straightforward approach is the following. First we truncate the 
Hamiltonian (5) at a finite polynomial order in the actions. This is legitimate if the 
radius g of the domain is small, due to the well known properties of Taylor series. 
However, the Fourier expansion of every term Hg still contains infinitely many contri- 
butions. Here we take advantage of the exponential decay of the Fourier coefficients of 
analytic functions and of some algebraic properties of the Poisson brackets. Precisely, let 
/(?) — /fc^*^'^''^^' dependence of the coefficients fk on the actions is unrelevant 
here. Then the exponential decay of the coefficients means that l/fcl < Ce-I'^l'^ with 
some positive constants C and a. Thus, having fixed a positive integer K we truncate 
the Fourier expansion as f(q) = 'Yli\k\<K fk^^^^''^^ remove all Foiuier modes 
\k\ > K. This is allowed because the exponential decay assures that the neglected part 
is small. The interested reader may find a more detailed discussion about this method 
of splitting the Hamiltonian in [18]. 

Coming back to our problem, we include in Hg (p, q) all Fourier coefficients with 
\k\ < sK, so that Hs{p, q) is a homogeneous polynomial of degree s + 1 in the actions 
p and a trigonometric polynomial of degree sK in the angles q. The algebraic property 
mentioned above is that such a splitting of the Hamiltonian is preserved by the Lie series 
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algorithm that we apply through all our calculations. This in view of the elementary fact 
that the Poisson bracket between two functions, fr and fs say, which are homogeneous 
polynomial of degree r + 1 and s + 1, respectively, in p and trigonometric polynomials 
of degree rK and sK, respectively, in q produces a new function of degree r + s + 1 in 
p and (r + s)K in q. 

A final remark concerns the estimate of the remainder J> in (8), which is an infinite 
series, too. Here we just calculate the first term of the remainder, namely the term of 
degree r + 1, and multiply its norm by a factor 2. This factor is justified in view of the 
fact that the analytical estimates of the same quantities involve a sum of a geometric 
series which, for g small enough, decreases with a ratio less than 1/2. Here a natural 
objection could be that for some strange reason the norm of the remainder at some 
finite order could be smaller than predicted by the analytical estimates. However, it is 
a common experience that after a few perturbation steps the norms of the functions 
take a rather regular behaviour consistent with the geometric decrease predicted by the 
theory. Thus, our choice appears to be justified by experience. 

As a final remark we note that our way of dealing with the truncation is the most 
straightforward one, but it is not the sole possible. Other more refined criteria may 
be invented, of course, which may take into account the most important contribution 
while substantially reducing the number of coefficients to be calculated. In this sense 
our direct approach should be considered as a first attempt to check if the concept of 
Nekhoroshev 's stability may be expected to apply to our Solar System. Although being 
unable to produce rigorous results in a strict mathematical sense, we believe that our 
method gives interesting results in the spirit of classical perturbation methods. 

4. Application to the planetary problem of three bodies 

Applying the theories of Kolmokorov and Nekhoroshev to the planetary problem is not 
straightforward, due to the degeneration of the Keplerian motion. In order to remove 
such a degeneration, a lenghty procedure is needed; this essentially requires a suitable 
adaptation of the canonical coordinates, paying a very particular care to the secular 
ones (to appreciate some deep point of view about this problem, see, e.g., [2] and [35]). 

In our approach the difl&culty shows up in the part concerning the application of 
Kolmogorov's theory. Once a Kolmogorov torus has been constructed, then there is no 
extra difficulty in applying the method of sect. 3, due to the fact that the method is 
local. In the present section we give a brief sketch of the procedure for the construction 
of a Kolmogorov torus. The complete procedure is described in [27] and [28], to which 
we refer for details. 

Following a traditional approach, we first reduce the integrals of motion (i.e. the 
linear and angular momenta); therefore, we separate the fast variables (essentially the 
semimajor axes and the mean anomalies) from the slow ones (the eccentricities and the 
inclinations with the conjugated longitudes of the perihelia and of the nodes). This is 
usually done in Poincare variables by writing a reduced Hamiltonian of the form 

(10) if«(A, A, ry) = (A) + fiF^'^ (A, A, ^, ry) 
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with 

fx = max{mi / niQ , m2 / mo } 

where mo is the mass of the star, nii and m2 are the masses of the planets, Aj, Xj are 
the fast variables and ^j, rjj are the slow (Cartesian- like) variables. Here, obviously, the 
values of the index j = 1 , 2 correspond to the internal planet and to the external one, 
respectively. 

On this Hamiltonian we perform a procedure which is the natural extension of the 
one devised by Lagrange and Laplace in order to calculate the secular motion of the 
eccentricities and the inclinations and the conjugated angles. 

The first step is the identification of a good unperturbed invariant torus for the fast 
angles A, setting for a moment the slow variables ^, rj to zero. Here is a short description, 
(i) Having fixed a frequency vector n* G R^, we determine the corresponding action 
values A* corresponding to a torus which is invariant for an integrable approxima- 
tion of the system, where the dependency on both the fast angles A and on the 
secular coordinates ^, r] is dropped. This can be done by solving the equation 



d{H 



R\ 



aA, 



A=A* 



= <, j = l,2. 



Here {H^)x = d\\d\2 is the average of the Hamiltonian with 

respect to the fast angles. The explicit value of n* is chosen so that it refiects 
the true mean motion frequencies of the planets (see next section for our values). 
Having solved the previous equation with respect to the unknown vector A*, we 
expand in power series of A — A* . With a little abuse of notation we denote 
again by A the new variables. 

(ii) We perform two further canonical transformations which make the torus A = ^ = 
?7 = to be invariant up to order 2 in the masses. Indeed, these changes of coordi- 
nates are borrowed from the Kolmogorov's normalization algorithm, but we look for 
a Kolmogorov's normal form with respect to the fast variables only, considering the 
slow ones essentially as parameters, although they are changed too. More precisely, 
we determine generating functions of the form Xj(A,A, ^, 77) = A.^~^gj{\^,ri) for 
J = 1,2 , where gj{X,C,f]) includes a finite order expansion both in Fourier modes 
with respect to the fast angles A and in polynomial terms of the slow variables ^,77. 
The aim of this step is to reduce the size of terms independent of or linear in the 
fast actions so that it is of the same order as the rest of the perturbation. We denote 
by the resulting Hamiltonian, which is still trigonometric in the fast angles A 
and polynomial in A, ^,77. 

The next goal is to determine a good invariant torus for the slow variables ^,r] . To this 
end we combine the classical Lagrange's calculation of the secular frequencies with a 
Birkhoff's procedure that takes into account the nonlinearity. 

(iii) We consider the secular system, namely the average {H'^) of the Hamiltonian 
resulting from the step (ii) above. Acting only on the quadratic part of the Taylor 
expansion of (H^) in ^ , 7] we determine a first approximation of the secular fre- 
quencies, and transform the Hamiltonian so that its quadratic part has a diagonal 
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Table 1. Physical parameters for the Sun-Jupiter-Saturn system taken from JPL 
at the Juhan Date 2451220.5 . 





Jupiter (j = 1) 


Saturn (j = 2) 


mass rrij 


(27r)Vl047.355 


(27r)V3498.5 


semi-major axis aj 


5.20092253448245 


9.55716977296997 


mean anomaly Mj 


6.14053316064644 


5.37386251998842 


eccentricity Cj 


0.04814707261917873 


0.05381979488308911 


perihelion argument loj 


1.18977636117073 


5.65165124779163 


inclination ij 


0.006301433258242599 


0.01552738031933247 


longitude of the node Qj 


3.51164756250381 


0.370054908914043 



form. This part of the calculation follows the lines of Lagrange's theory, but the 
calculation is worked out at the second order approximation in the masses. The 
diagonalization of the quadratic part requires a linear canonical transformation, 
which is a standard matter. Thus the quadratic part in ^, rj of the resulting Hamil- 
tonian has the form | + ^f)' where u are the secular frequencies and we 

denote again by ^, the slow variables. 

(iv) We perform a Birkhoff 's normalization up to order 6 in ^, i]. This gives a normalized 
secular Hamiltonian which in action-angle variables — ^'ilj cos ipj , r]j = 
■\/'2.Ij sin ipj takes the form 

H"" = vl + (/) + (/) + F(A, /, ^) , 

where h'^^^ and h^^^ are polynomials of degree 2 and 3 in /, respectively. This 
step removes the degeneration of the secular motion, thus allowing us to take into 
account the nonlinearity of the secular part of the problem. 

(v) Having fixed the slow frequencies g* so that they reflect the true frequencies of the 
system, we determine a secular torus /* corresponding to these frequencies, by using 
the integrable approximation of . This is done by solving for / the equation 

The values A* and /* so determined provide the first approximation of the Kolmogorov's 
invariant torus. Reintroducing the fast angles and performing on the original Hamilto- 
nian all the transformations that we have done throughout our procedure (i)-(v) 
wc get a Hamiltonian of the form (2) which is the starting point for Kolmogorov's nor- 
malization algorithm. After a number of Kolmogorov's steps the Hamiltonian takes the 
form (4), thus giving a good approximation of an invariant torus with frequencies n* 
and g*. The latter form is precisely the output of the calculation illustrated in [17], and 
by removing all terms which are independent of or linear in the actions p it provides a 
Hamiltonian as that in (5). This is the starting point for our algorithm evaluating the 
stability time in the neighbourhood of the invariant torus. 
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Table 2. The frequencies of the unperturbed torus in the SJS system correspond- 
ing to the initial data and physical parameters in table 1. The values are calculated 
via frequency analysis on the orbits obtained by direct integration of the equations 

for the yiroblcni of tlirco bodies. 





Jupiter 


Saturn 


fast frequencies 


nl = 0.52989041594442 


n* = 0.21345444291052 


secular frequencies 


gl = -0.00014577520419 


g* = -0.00026201915143 



5. Application to the Sun— Jupiter— Saturn system 

We come now to the application of our procedure to the SJS system. Let us first de- 
fine the model. We consider the general problem of three bodies with the Newtonian 
potential. Thus, the contribution due to the other planets of the Solar System is not 
taken into account in our approximation. The expansion of the Hamiltonian is a classi- 
cal matter, so we skip the details, just recalling that all the expansions have been done 
via algebraic manipulation, using a package developed on purpose by the authors. 

The choice of the model plays a crucial role in determining the frequencies of the 
torus, that we calculate by integrating the Newton equations for the problem of three 
bodies and applying the frequency analysis (see, e.g., [22]) to the computed orbit. As 
initial data we take the orbital elements of Jupiter and Saturn as given by JPL-*^ for the 
Julian Date 2451220.5 . This is the point where the connection with the physical param- 
eters of our Solar System is made. The physical parameters and the orbital elements 
are reported in table 1. The calculated frequencies are given in table 2. 

The choice of the Julian Date 2451220.5 in order to set the initial data is completely 
arbitrary, of course, its sole justification being that such data are directly available from 
JPL. Choosing different dates or different determinations of the planet's elements could 
lead to a slightly different determination of the frequencies, and so also of the invariant 
torus. However, we emphasize that the aim of the present work is precisely to give a long 
time stability result which applies to a neighbourhood of the invariant torus. The size 
of such a neighbourhood should be large enough to cover the unavoidable uncertainty 
in determining the initial data for the SJS system. This is a delicate matter, of course, 
because the JPL data reflect the dynamics of the full Solar System, while our study is 
concerned only with the model of three bodies. However, we may get some hint on the 
size of the uncertainty precisely by looking at the JPL data. 

As everybody knows, the initial positions and velocities of the planets are usually 
taken from the Development Ephemeris of the Jet Propulsion Laboratory (for short, 
JPL's DE). There are several sets of these ephemerides, each version of them being 
based on more and more observational data, which take benefit from the improvement 
of the techniques. Thus, each new version of the JPL's DE is expected to improve the 
precision of the data with respect to the older ones and, then, one can approximately 

^ The data about the planetary motions provided by the Jet Propulsion Laboratory are 
publicly available starting from the webpage http://www.jpl.nasa.gov/ 
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Table 3. Estimates of the uncertainties on the initial values of the canonical coor- 
dinates (A, A, ^, r/) . These evaluations are derived from the comparison of different 
sets of JPL's DE. 





AA, 


AA, 






Jupiter (j = 1) 


1.8 X 10-6 


6.6 X 10-^ 


1.1 X 10-^ 


2.8 X 10-6 


Saturn {j = 2) 


1.7 X 10-6 


3.0 X 10-^ 


3.3 X 10-6 


3.2 X 10-6 



Table 4. Maximal discrepancies about the orbital elements of the SJS system 
between a numerical integration and the semi-analytic one, that is based on the 
construction of the invariant torus corresponding to the frequencies values given 
in table 2. The maximal relative errors on the semi-major axis aj and on the 
eccentricities ej are reported here for both Jupiter (corresponding to j — 1) and 
Saturn (i.e., j = 2); the same is made also for the maximal absolute errors on the 
"fast angle" Xj = Mj + Uj and on the perihelion argument Uj . In the present case, 
the comparisons are made starting from the initial conditions given in table 1 and 
for a time span of 100 Myr. 







Maxi{|AA,(t)|} 




MaXi{|Aa;,-(t)|} 


Jupiter 


1.5 X 10-6 


5.0 X 10-4 


1.3 X 10-3 


1.3 X 10-3 


Saturn 


6.8 X 10-6 


1.1 X 10-^ 


4.3 X 10-3 


7.3 X 10-3 



evaluate the error of the older versions by comparison with the most recent one (see 
the initial discussion in [39]). The positions and velocities of the planets given by five 
different sets of JPL's DE are listed in table 15 of Standish's paper [38]. For each kind 
of these data we can determine a narrow interval containing all of them and we can 
calculate the Keplerian orbital elements corresponding to the extrema of such intervals. 
By applying all the necessary transformations we translate these data into uncertainties 
for the Poincare's canonical coordinates (A, A, ^, ij) that have been used in order to write 
the Hamiltonian (10). These uncertainties are reported in table 3. This provides us with 
a first approximation of the neighbourhood of our initial data that contains all JPL's DE 
reported in Standish's paper. We should now apply all the canonical transformations 
needed in order to construct an invariant torus close to the SJS orbit. However we remark 
that all such transformations are very smooth, being analytic, volume preserving, and 
most of them are close to identity, so that they add just a small correction with respect 
to the data in table 3. Thus we may confidently expect that at some time the phase 
space point representing the position of the SJS system lies in a neighbourhood of our 
approximated invariant torus the size of which is evaluated to be 0(10-6) for the fast 
actions and O{10~^) for the secular coordinates. 

Let us now come back to the actual calculation. The Kolmogorov's normal form has 
been computed up to order 17, with the generating function exhibiting a good geometric 
decay. Furthermore, we have compared the orbit on the approximate invariant torus with 
that produced by a direct numerical integration of the equations of motion, thus finding 
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a quite good agreement between them, as shown in table 4. Here, we omit the details 
about these lengthy calculations, since a complete report has been already given in [28]. 

The calculation of Kolmogorov's normal form produces a Hamiltonian which is an- 
alytic in the neighbourhood of the approximated invariant torus. Our program performs 
the calculation of this Hamiltonian with the polynomial series in the actions truncated 
at order 3 and the trigonometric series truncated at order 34 (see [28] for more details). 
On this Hamiltonian we would like to apply the procedure of sect. 3. However, a major 
obstacle raises up: the number of coefficients in the series that we have calculated is 
more than 7 100 000. Such a huge number of coefficients can not be handled in a Birkhoff 
normalization procedure. For, referring to the discussion at the end of sect 3 we should 
set the parameter K for the truncation of trigonometric series to 34, thus getting a 
truncation at trigonometric degree 68, 102, . . . , 34r, ... at successive order. A rough 
estimate of the number of generated coefficients shows that we shall soon run out of 
memory and of time on any available computer. Thus, we must introduce some further 
approximation. 

In view of the considerations above we decided, as a first approach, to strictly follow 
the truncation scheme illustrated at the end of sect. 3 by just lowering the value of K. 
We report the results of this ffist attempt, which in our opinion appear already to be 
interesting. Thus we expand the Hamiltonian in the form 

H{p, q) = {uj,p) + i?i {p, q) + H2{p, q) 

by keeping in Hi all terms of degree 2 in the actions p and K in the angles and in H2 
all terms of degree 3 in the actions p and 2K in the angles q. The Birkhoff normalization 
produces a Hamiltonian of the form 

H = {uj,p) + Zi{p) + . . . + Zr{p) + J>+i(p, q), 

where J>+i denotes the term of degree r + 2 in the actions p and (r + 1)K in the angles, 
i.e., the first term of the remainder. With a suitable choice of K, this considerably 
reduces the number of coefficients in the expansions thus enabling us to perform the 
calculation on a workstation. We emphasize however that the algorithm is a general one 
so that in principle it can be applied to the full Hamiltonian or, better to a Hamiltonian 
obtained by removing all coefficients which are very small and will likely not produce 
big coefficients (due to the action of small denominators) during the calculation of the 
Birkhoff's normal form. The rest of the calculation closely follows the discussion in 
section 3, so we come to illustrating the results. We performed the calculation with two 
different values of K, as given in the following table. 

K r # of coefficients 
4 5 2 494000 
6 4 3 380000 
This shows in particular the dramatic increase of the number of coefficient in the 
remainder J-r+i (third column), which imposes strong constraints on the choice of the 
normalization order r. 

A quite natural objection could be raised here. Since the most celebrated resonance 
of the SJS system (i.e., the mean motion resonance 5 : 2) has trigonometrical degree 7, 
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it seems that some of the main resonant terms arc neglected because of our choice of 
K . This is actuaUy not the case, due to a technical element that we have omitted in 
the previous section in order to make the discussion simpler. Our sequence of transfor- 
mations includes a unimodular linear transformation on the angles, and so also on the 
frequencies. The action on the frequencies changes the resonance 5 : 2 into a 3 : 1 one, 
which is of order 4. Thus, setting X > 4 as we did throughout all our calculations is 
enough in order to include the main resonant terms. The interested reader will find a 
detailed discussion of this point in sect. 3 of [28]. 

Let us now come to the results. In panel (a) of fig. 1 we report the results for K — 4. 
The crosses give the estimated stability time for the Birkhoff 's normal form at order 5; 
the dashed line gives the estimated time when the Birkhoff 's normal form is truncated 
at order 4, thus showing how relevant is the improvement when a single normalization 
order is added. 

We can now come back to the estimate about the escape time T = T{qq) . Looking 
at panel (a) of fig. 1, one can remark that we have an estimated stability time of 
IQi^^ years, that is approximately equal to the age of the universe, for a neighbourhood 
of initial conditions of radius 10~^ in actions. 

It may be noted that the stability curves exhibit a sharp change of slope around 
qq ~ lO"*^. This is because the optimal normalization order increases when the radius 
is decreased. Actually, further changes of slope should be expected for smaller values of 
£>o, but due to computational limits such changes can not appear in our figure, because 
the optimal order exceeds the actual order of our calculation. Thus, our estimate of the 
stability time should be considered as a very pessimistic lower bound. 

Moreover, the behaviour of the plots in fig. 1, clearly shows that the estimate of 
the escape time can be substantially improved if smaller values of the radius can 
be considered. Recalling that our estimate of the size of the neighbourhood in action 
variables is calculated from the discrepancy among different sets of JPL's DE data, we 
may aflfirm that our neighbourhood roughly covers such a width, which is tabulated in 
Standish's paper quoted above. 

If we try a better approximation of the Hamiltonian, setting = 6, then we are 
forced to stop the Birkhoff 's normalization at order 4, thus making the results definitely 
worse. The data for the estimated stability time are plotted in panel (b) of fig. 1. One 
sees that the estimate becomes comparable with the age of the universe only in a 
neighbourhood of initial conditions slightly larger than 10~^. However, if we compare 
the curve in panel (b) with the dashed curve in panel (a) we see that we shall likely get 
substantially better results if we could compute the normal form at order 5. 

Thus, our rough approximation gives results which apply to a set of initial data for 
the SJS system which is of the same order of magnitude as the uncertainty in JPL's 
data. We also emphasize that our evaluation of is based on observational data which 
are presently older than 25 years; we expect that this is quite pessimistic with respect 
to the features of more recent JPL's DE. 
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Figure 1. The estimated stability time, (a) results for K = 4 and BirkhofF normal- 
ization order 5 (crosses) and 4 (dashed curve), (b) results for = 6 and Birkhoff 
normalization order 4. 



6. Conclusions 



We have developed an effective method to compute the Kolmogorov's normal form for 
the problem of three bodies, and have successfully applied it to the SJS problem, using 
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the data of our Solar System. Next, we have shown that a calculation of the stability 
in Nekhoroshev's sense of orbits with initial point in the neighbourhood of the torus is 
possible, at least if one accepts to make some strong truncation in the expansion of the 
Hamiltonian. A rather strong truncation allows us to get an estimated stability time 
comparable with the age of the universe in a neighbourhood of the torus that will likely 
contain the actual initial data of the SJS system. 

The natural question is whether such results will remain valid if we add more and 
more terms in the Hamiltonian, thus making our approximation better. Answering such 
a question is presently beyond our limits, but in our opinion deserves to be investigated. 
Some improvement may be obtained by using more computer power, e.g., by performing 
our calculation on a cluster of computers. However, substantial improvements require 
also a refinement of our analytical techniques in order to be able to evaluate the error 
induced by our truncations, thus allowing us to introduce better computation schemes 
and to evaluate the reliability of our approximations. 
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